Proceedings of the 11th International Conference 
on Computational and Mathematical Methods 
in Science and Engineering, CMMSE 2011 
26-30 June 2011. 



An 0{N^) implementation of Hedin's GW approximation 



Peter Koval^, Dietrich Foerster^ and Daniel Sanchez-Portal^ 

1 Centra de Fisica de Materiales CFM-MPC, Centro Mixto CSIC-UPV/EHU and 
DIPC, E-20018 San Sebastian, Spain 

^ CPMOH/LOMA, University of Bordeaux, France 

emails: koval.peter@gmail.com, d.foerster@cpmoh.u-bordeauxl.fr, 

sqbsapodSsq . ehu . es 



Abstract 



Organic electronics is a rapidly developing technology. Typically, the molecules 
involved in organic electronics are made up of hundreds of atoms, prohibiting a the- 
oretical description by wavefunction-based ab-initio methods. Density-functional 
and Green's function type of methods scale less steeply with the number of atoms. 
Therefore, they provide a suitable framework for the theory of such large systems. 

In this contribution, we describe an implementation, for molecules, of Hedin's 
GW approximation. The latter is the lowest order solution of a set of coupled 
integral equations for electronic Green's and vertex functions that was found by 
Lars Hedin half a century ago. 

Our implementation of Hedin's GW approximation has two distinctive features: 
i) it uses sets of localized functions to describe the spatial dependence of correlation 
functions, and ii) it uses spectral functions to treat their frequency dependence. 
Using these features, we were able to achieve a favorable computational complexity 
of this approximation. In our implementation, the number of operations grows as 

with the number of atoms N . 

Key words: Hedin's GW approximation, basis of dominant products, large 
molecules. 



1 Introduction 

The promising field of organic electronics deals with large molecules of several tens or 
even hundreds of atoms jlj. For instance, fullerene Cgo is a frequently used subunit in 
organic electronics and it alone consist of 60 atoms (see figure [T]) . 
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Each individual molecule may be used in a device in 
many different ways and there is an astronomically large 
number of different promising molecules. As in many cases 
there is a limited knowledge of the relevant physical pa- 
rameters, and it might be also interesting to explore the 
potential of candidate molecules theoretically, before these 
molecules has been actually synthesised. 

The geometry of large organic molecules can be re- 
liably predicted by density-functional theory (DFT)[3]. 
However, the properties of their excited states such as 
the energy of the highest occupied (HOMO) and lowest 
unoccupied molecular orbitals (LUMO), corresponding to 
adding and subtracting one electron from the system re- 
spectively, require a description of electronic correlations 
better than that provided by current functionals of DFT 
and its time-dependent counterpart, TDDFT. 

Such effects can be efficiently incorporated with the help of Hedin's method that is 
based on Green's function. Hedin's GW approximation for one-electron Green's func- 
tion is computationally cheaper than wavefunction-based methods, although it remains 
computationally more expensive than DFT and TDDFT within linear response. 

The goal of our work is to develop a practical algorithm for Hedin's GW approx- 
imation which is suitable for large organic molecules, allowing to access the excited 
states of such molecules. 



Figure 1: Ball and stick 
model of fullerene Ceo 
produced with XCrysDen 
package [2]. 



2 Theoretical framework for 

Electronic Green's function (propagators) are 
useful in condensed matter physics because 
many simple observables can be computed 
in terms of them. At the same time, such 
Green's functions remain simpler than many- 
body wavefunction. 

Hedin's GW is a useful approximation for 
the so-called self-energy S(r, r, w) that enters 
Dyson's equation for an interacting electronic 
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Figure 2: Feynman diagram of Dyson 
equation ([l]). 

propagator G{r, r' , oj) 



G-\r,r',u) = G^\ry,u:) - S(r,r,a;). (1) 

Here, the inversions must be understood in operator sense J (r, r", uj)G{r" , r', uj)dr" = 
5{r — r') and Go{r,r',uj) stands for Green's function where electron-electron interac- 
tions have been switched off. It is obtained from an effective one-particle Hamiltonian 



(w - H{r))Go{r,r',uj) = 6{r - r'). 



(2) 
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In this work we use a Kohn-Sham Hamiltonian [3], although Hartree-Fock Hamil- 
tonian also proves to be useful at this point [Ij. Hedin's GW approximation for the 
self-energy S(r,r,a;) reads 

S(r, r',t) = iGo(r, r' , t)Wo{r, r' , t). (3) 



It involves the non interacting electronic Green's 

function Go{r,r',t) and a screened Coulomb inter- ^Vq 
action TVo(^5 t). This approximation is a solution 
of a truncated version Hedin's equations [H [6]. The 
name of this approximation is taken from the simple 



form of the electronic self-energy S = iGW. t:^- o ^■ c 

l^igure S: I^eynman diagram of 

The screened Coulomb interaction Wq can be 



easily calculated in frequency domain using the so- 
called RPA approximation 



self-energy ([s]). 



Woir,r',oj)=[5{r-r"')-v{r,r")xo{r",r"',u)] v{r"',r'), (4) 

where v{r,r') = \r — r'\~^ is the bare Coulomb interaction. Here and in the following 
we assume integration over repeated spatial coordinates (r" and r'" in equation Q) 
on the right hand side of an equation if they do not appear on its left hand side. The 
screened interaction Q is the sum of the bare Coulomb interaction created by a point 
charge at r' , plus a correction due to the redistribution of charge induced in response 
to the total field [316]. The non-interacting response function xo('")'"';i) is related to 
the non-interacting Green's function 

iXo(r, r', t) = 2Go(r, r' , t)Go(r', r, -t), (5) 

where a factor 2 arises because of the summation over spin variable. 

As we mentioned already, 
we construct the non-interacting y y /^''N W 

Green's function using an ef- iVmrV^ = i.yv\/\xv2 + l Yirir!rTT*2 
fective Kohn-Sham Hamilto- 

Figure 4: Feynman diagram of screened Coulomb in- 
Hks = -Iv' + Vks, (6) action (§. 

6E 

Vks = Kxt + ^Hartroe + ^xc, where V^c{r) = T-p^- 

dn{r) 

Exc is a functional of the electronic density that includes the effects of exchange and 
correlation in an effective way. Its functional derivative Vxcir) is the so-called exchange- 
correlation potential and it must be subtracted from Ti(r,r',t) to avoid including the 
exchange-correlation interaction twice in equation ([s]). This is accomplished with the 
substitution 

$](r, r', t) ^ S(r, r', t) - 6{r - r')5(t)Kcc(r-) 
in Dyson's equation ([T]). 
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3 A basis set of localized functions 

Having the equations ( 1|3|4| 5) at hand we introduce a basis set of localized functions 
and rewrite the system of equations in the basis. We start with linear combinations of 
atom orbitals (LCAO) to represent the non-interacting Green's function GQ{r,r',t) 

Go{r,r',t) = Y,GUt)nr)f\r'), (7) 

ab 

where /""{r) are atom centered orbitals. The frequency (and time) dependence has been 
factorized in the last equation. The treatment of the frequency (and time) dependence 
by spectral functions will be explained in section [4| Inserting equation ([T]) into the 
equation we obtain 

iXo(r, r', t) = 2Y, G'i,{t)G',,{-t) nr)f{r) f\r')r{r'). (8) 

abed 

Products of localized orbitals such as f°'{r)f'^{r) appear in the last equation. Although 
a product of localized orbitals is also a localized function, such products do not form 
a suitable basis because they contain many collinear functions. Several methods have 
been proposed to construct more efficient basis to span the products of localized orbitals 
[HEllin]. Here we use a basis of dominant products mj that is constructed individually 
for each atom pair. The dominant products are identified as certain linear combinations 
of the original orbital products and they are free of any coUinearity within a given 
atom pair (with respect to a given metric, here we have used the Coulomb metric). 
Moreover, the original orbital products can be expressed as linear combinations of 
dominant products 

nr)f\r) = V^'F^{r). (9) 

The three-index coefficient will be referred to as the product vertex. The product 
vertex is local or sparse by construction and indeed the locality of our construction is 
its main characteristic. 

Considering Dyson's equation ([T]), we arrive at its tensor counterpart 

Gabiuj) = G°,(^) + Gaa'{uj)^^'''\i0)GU^), (10) 

where matrix elements of the self-energy S°^(a;) must be used 

^^\uj) = jj r{r)^{r, r', uj)f\r') d^rSr' . (11) 



Calculating the matrix elements of the self-energy by equation ([s]) and using ([T]) 
for the non interacting Green's function, we arrive at 

T,-\u:) = i^Gl,(t) / nr)r' {r)W,{r,r\t)f [r')f\r')d\d\' . (12) 

a'b' 

Using the identity ([9]), we rewrite the latter equation as 

^^\u:) = \G\,{t)VfW^''{t)V^'\ (13) 
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where the matrix elements of the screened Coulomb interaction appear 

W^''{t) = jjF^'{r)Wo{r,r\t)F^{r')(frd^r'. (14) 

Finally, the equation Q gives rise to the corresponding tensor expression 

Wi;^{u) = (5^, - v^^'x%{^))-\'''''. (15) 

The last expression can be elucidated by developing the operator [1 — fXo]~^ in a 
geometric series [1 — vxo]~^ = 1 + vxo + ^^Xo'^Xo + vxqvxqvxo ■ ■ ■ The expressions 



10[ ), (13) and (15) are tensor counterparts of Hedin's equations in coordinate space 
(jl), ([3]) and ([4]), correspondingly. In the next section, we will present our method 
for treating the frequency (and time) dependence of these tensor equations. 

4 Spectral function technique 

Because of the discontinuities of the electronic Green's functions, a direct, straightfor- 
ward and accurate computation of the response function ([s]) is practically impossible 
both in the time domain and in the frequency domain. However, one can use an 
imaginary time technique [12] or spectral function representations to recover a com- 
putationally feasible approach. In this work, we use spectral function techniques and 
rewrite the time ordered operators as follows 

Jo J -00 

POO pO 

xlAt) = -m / d5a+,(s)e--* + i^(-t) / dsa'^Me-'^'; 

Jo J~oo 
POO f'O ^ ' 

W^'^it) = -Wit) ds-t1^{s)e-'''* + \e{-t) ds-i't%s)e-'''- 
Jo J -00 

/•oo fO 

T.''\t) = -Wit) dsal\s)e-''' + ie{-t) dsal\s)e-''\ 
Jo J -00 

where "positive" and "negative" spectral functions define the whole spectral function by 
means of Heaviside functions 0[t). For instance, the spectral function of the electronic 
Green's function reads 



Transforming the first of equations ( 16 ) to the frequency domain, we obtain the familiar 



expression for the spectral representation of a Green's function 



00 



Here e is a small line-broadening constant. In practice, the choice of e is related to the 
spectral resolution Aw of the numerical treatment. 
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As first application of representations ( |16[ ), we derive the spectral function of the 
non interacting response a^y{s) using equation ([5]) as a starting point 



= jj Wc(^l)K>^a(-«2)5(si + S2- s)ds,dS2. (18) 

Here, the convolution can be computed with fast Fourier methods and the (time- 
ordered) response function x^ui^) be obtained with a Cauchy transformation 



oo 



xlA^) = X^i.(-w) + x;i'^(w), where xM^') = I ~rh ' ^^^^ 

Jo u! -\-ie — s 

The calculation of the screened interaction Wq'^{uj) must be done with functions, 



rather than with spectral functions, because of the inversion in equation (15). The 
spectral function of the screened interaction ^^^[uj) can be easily recovered from the 
screened interaction itself [^. Since Im^-q^^— ^ is a representation of Dirac (5-function 

when e goes to zero, then j^'^{uj) = ImWo''^ (uj). 

Deriving the spectral function (t{lo) of the self-energy, we arrive at 

af{s)= / / 6{si + S2-s)V;:-'p^,,,{si)V,'''^1''{s2)dsids2, (20) 
Jo Jo 

rO rO 



I I 5{s, + S2- s)V^''' p-,,,{s,)V,''V-{s2)ds,ds2. 

J —oo J —oo 



These expressions show that the spectral function of a convolution is given by a con- 
volution of the corresponding spectral functions. 

4.1 Discretization of frequency-dependent quantities 



The spectral functions in equation (18) are merely a set of poles at (eigen) energies E 



E>0 E<0 

Here the eigenvectors diagonalize the corresponding Kohn-Sham Hamiltonian 

Tjab-w-E TpaabvE 

where the Hamiltonian and the overlap matrices of atomic orbitals f'^{r) enter 

jjab ^ j f^(^r)H{r)f\r)d^r, and S''^ = J r{r)f\r)d^r. (22) 

In practice, we use the SIESTA package jl3j that gives the orbitals /"(r), eigenvectors 
and eigenvalues E for a given molecule as the output of a DFT calculation. 



The use of fast Fourier techniques for convolution, for instance in equation (18), 



requires that the spectral functions />^(cj), p^^i^) known at equidistant grid points 
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Uj = jAuj,j = —Ni^ . . . Ni^, rather than at a set of energies resulting from a diagonal- 
ization procedure. The solution for this problem (discretization of spike-like functions) 
is known and well tested [T3] . We define a grid of points that covers the whole range of 
eigen energies E. Going through the poles E, we assign their spectral weight X^X^ 
to the neighboring grid points n and n + 1 such that Un < E < uJn+i according to 

the distance between the pole and the grid points Pn,ab = " X^X^ , Pn+i,ab = 

1 —Pn,ab- Such a discretization keeps both the spectral weight and the center of mass 
of a pole. It also reduces the number of operations that are needed to calculate the non 
interacting response function x^^ui^)- This is so because the number of frequencies A'^^ 
can be kept small (typically a few hundred points) even for large molecules while the 
number of states A'ci-b grows linearly with the size of the system. 

4.2 The second window technique 

The discretization of spectral weight helps to control the computational complexity 
for large molecules. However, we are actually interested in the properties of low lying 
levels (HOMO and LUMO and several levels below and above). At first sight one might 
think that one could neglect the contributions of high energy spectral weights in the 
Cauchy transformation. However, neglecting the high energy spectral weight actually 
results in a wrong real part of the functions. Fortunately, the high energy spectral 
weight tolerates a coarser grid |8j. Therefore, we calculate each spectral function twice: 
once with a higher resolution in a low frequency range, and a second time with a lower 
resolution but in the whole range. The Cauchy transformation for such a two-window 
representation must be modified as follows 



X^^(a; + iesmall) = / ds — + / + ] ds 



, small window/ , , ■ \ , 

X,j.u + lesmallj + 



X iO -\- lEsmall — S \J-A J\ J + l^large — S 

large window/, , i ;^ 
tiu + leiarge 



, large window/, , i ;^ \ /'OQ^ 
Xnu + leiargej , , „ . • (^-^j 



truncated spectral function 



After the calculation of spectral functions in both windows, we truncate the spectral 
function in the second window in the range ... A, do Cauchy transform of both spectral 
functions and update (by a linear interpolating procedure) the function in the first 
window with the truncated function from the second window. 

We use the second window technique both for the non interacting response function 
X^j,(w) and for the self-energy S"^(a;). 

5 Non-local compression of the product basis 

The basis of dominant products is optimal within a given atom pair, but unfortunately, 
there is still a lot of collinearity between dominant products belonging to different 
pairs. This collinearity is an indication that the size of the product basis can be 



strongly reduced. Even for the molecules of modest size considered in Section 7.1 the 
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basis set of dominant product becomes so large that hampers the storage of the (non- 
interacting) response function (19) and slows down the inversion in the calculation of 
the screened interaction (15). In order to improve the situation we perform a non- 
local (global) contraction of the basis of dominant product. We start by considering 
a sum-over-states expression for the non-interacting response function in the basis of 
dominant products 



XU-) = where V^^^ = Xf <Xf . (24) 



The response x^ui^) built up from vectors that represent electron-hole pair 

excitations. One can use these vectors to identify important directions in the space 
of dominant products. The number of electron-hole pairs EF grows as with the 
molecular size while the dimension of dominant product basis is 0{N) by construction 
(due to the localization of the basis orbitals). Therefore, one has to limit the set 
of electron-hole pairs EF from the beginning to keep the efficiency of the algorithm, 
particularly if one uses a diagonalization-based procedure for generating the (globally) 
optimal basis. Because of the inherent limitations of LCAO to represent high energy 
features, and the fact that we are mainly interested in the lowest lying excitations, we 
choose 0{N) low-energy electron- hole pairs 

{X;} = subset of {V^^} limited by \E - F\ < ^threshold, n = l... iV^^nk- (25) 

After the initial selection according to the energy criterion \E — F\ < i^thrcshoidi we 
define a metric g'"" 



9"'" 



whe.e.-.//F'.M|.-.'r'F-(/).W.'. (26) 



After diagonalizing the metric ^'""Cn = -^^mi can identify important directions (like 
in the construction of the basis of dominant products [HI [T^] ) by building linear combi- 
nations of the original vectors and by choosing only eigenvectors with eigenvalues 
above a suitable threshold value 



^ K^i/^- (27) 

These linear combinations can be used to expand the original response function x^^y{oj) 
in terms of fewer functions 

X%{^) = Z-;xln{^)Z:. (28) 
In order to express x^ni^) terms of x^^ui!^) multiply equation (28) with Z'^v^'^ 



from both sides and notice that Z'^v^^ Z'^ = Z'^Zn = 5™. Therefore, the response 
function can be "compressed" by using basis vectors Z^ = v^'^Z^: 

Xlni^) = Z^x%{^)K- (29) 
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The particular choice of the Coulomb metric v^'^ in equation (26) simplifies the 
computation of the Coulomb screened interaction (15). We can rewrite the equation 
(15) in terms of a Taylor series 

Wr = V^" + v'^'^'x^v"'" + v'^'^'x%v'''>^"x%.v''"'' + ■■■ (30) 



Inserting here the response function x^u according to equation (|28|) and recalling the 
identity ZJ^ZH = 5™, one arrives at 



= V'^'' + Z^xL [Srn + Xrn + XrsXL + " " " ] 



(31) 



"T ^mAmn ^n! "'^^^^^ Amn — \^mk Xmk) A.kn- 



At this point it should be also be noted that the self-energy S^^(a;) that corresponds 
to the instantaneous part of the screened interaction v^^ is computed separately (Gj E] 
without any non local compression. 



6 Computational complexity of the algorithm 

The number of mathematical operations spent in different parts of the approach pre- 
sented above can be estimated if the dimensions of the corresponding matrices are 
known. The numbers that determine the complexity of the algorithm are the number 
of atomic orbitals A'orb! the number of dominant functions A'prod and the number of 
frequencies N,^. The number of orbitals and the number of dominant products are 
proportional to the number of atoms N by construction. The number of frequencies af- 
fects the run time linearly, but it is independent of the number of atoms. The non-local 
basis of section pi can be constructed in 0{N^) operations because A^rank in equation 



(25) can be kept proportional to number of orbitals. In practical calculations we have 
found that converged results are achieved with Arank ~ 5Aorb- For large molecules, 
the number of important eigenvectors Agubrank after dropping small eigenvalues A in 



equation (27) is approximately Aorb- No part of the algorithm scales worse than O(A^) 
[8]. There are several portions of the code where 0{N^) operations are needed. How- 
ever, only two of them have an appreciable impact on the run time: the computation 
of the response function and the computation of the self-energy. Both of them scale as 
^(-^prod-^subrankAtj) and give rise to an overall O(A^) scaling of the run time. 



7 Applications to organic molecules 

The methods described in the previous sections were carefully tested on several molecules. 
In this paper, we present two examples: calculations of HOMO and LUMO levels of 
three aromatic hydrocarbons (benzene, naphthalene and anthracene) and a calculation 
of the HOMO and LUMO levels of fullerene Ceo- 
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7.1 Aromatic hydrocarbons 



From the interacting Green's function Gab{^) we calculate the density of states (DOS) 
p{uj) = —S"'^lm.Gab{^)/'^ and we then determine the positions of the HOMO and 
LUMO levels from the DOS. 

The results of this procedure for aro- 
matic hydrocarbons are collected in table 
[!} One can see that our (LDA-l-GoWo) 
approach delivers qualitatively correct pre- 
dictions for the ionization potentials (IP) 
and electron affinities (EA) of benzene and 
naphthalene (donors) and anthracene (ac- 
ceptor). On the other hand, the LUMO 
of the underlying DFT calculation is always 
below the vacuum level. The calculations 
have been done on top of DFT-SIESTA cal- 
culations. We used pseudo potentials of 
Troullier-Martins type [17] and the Perdew- 
Zunger exchange-correlation functional [18j . 
We found that rather extended atomic or- 
bitals must be used to achieve converged re- 
sults in our GW approach. The energy shift 
parameter [19], that controls the spatial ex- 
tension of atomic orbitals has been set to 3 meV for benzene, and to 20 meV for naph- 
thalene and anthracene. The spectral functions have been discretized in two energy 
windows, with each window containing N^^ = 64 frequency points. 



Picture 


IP, eV 


EA, eV 




8.82 


-1.43 




(9.25) 


(-1.12) 




7.58 


-0.15 




(8.14) 


(-0.19) 




6.87 


0.73 




(7.44) 


(0.530) 



Table 1: The ionization potentials (IP) 
and electron affinities (EA) of benzene, 
naphthalene and anthracene. Experi- 
mental data [16] are given in brackets. 



7.2 Fullerene C 



Source 


IP, eV 


EA, eV 


Our LDA+GoWo 
Experimental p6j 


7.33 
7.58 


2.97 
2.65 



Table 2: The IP and EA of fullerene 
Geo calculated with our method and 
corresponding experimental data. 
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The fullerene Ggo and its derivatives are very 
popular ingredients in organic semiconductors 
and extensive experimental data and theoretical 
computations are available for the basic fullerene. 
We found that our LDA-|-GoVFo results are in 
very good agreement with experimental data (see 
table [2|. The computational parameters of this 
calculation are the same as in subsection |7.1[ while the energy shift parameter is cho- 
sen to be 15 meV. The number of frequency points was chosen rather large A'^^ = 128 
and the calculation has been done with 8 cores of a Nehalem machine (Intel(S)E5520 
2.27GHz, Gache 8M/DDR3 RAM 24GB). The current version of the code consumed 
26 hours of wall clock time. 

A comparison of DOS calculated with DFT LDA Hamiltonian and with our LDA-|-Go Wq 
approach is shown in figure [5j Such a result is a typical when Hedin's GW approach is 
applied on top of a LDA calculation. GW HOMO has lower energy than DFT HOMO. 
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Therefore, the change density n(r) will be more localized in the GW calculation. GW 
LUMO has higher energy than DFT LUMO. Therefore, the change density n(r) will 
be more delocalized in the GW calculation. 

8 Conclusions 

We have described our approach to Hedin's 
GW approximation for finite systems. This 
approach allows to compute the interact- 
ing Green's function on a frequency grid. 
The density of states is our output and it 
provides HOMO and LUMO levels in rea- 
sonable agreement with experiment. The 
complexity of the approach scales with the 
third power of the number of atoms, while 
the needed memory scales with the second 
power of the number of atoms. These fea- 
tures make our approach suitable for treat- 
ing the large molecules that are used in or- 
ganic semiconductors. 
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